Quantum phase transitions in the honeycomb-lattice Hubbard model 
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Quantum phase transitions in the Hubbard model on the honeycomb lattice are investigated in the 
variational cluster approximation. The critical interaction for the paramagnetic to antiferromagnetic 
phase transition is found to be in remarkable agreement with a recent large-scale quantum Monte 
Carlo simulation. Calculated staggered magnetization increases continuously with U and thus we 
find the phase transition is of a second order. We also find that the semimetal-insulator transition 
occurs at infinitesimally small interaction and thus a paramagnetic insulating state appears in a wide 
interaction range. A crossover behavior of electrons from itinerant to localized character found in 
the calculated single-particle excitation spectra and short-range spin correlation functions indicates 
that an effective spin model for the paramagnetic insulating phase is far from a simple Heisenberg 
model with a nearest-neighbor exchange interaction. 

PACS numbers: 71.10.Fd, 71.27.+a, 71.30.+h 



I. INTRODUCTION 

The correlation induced metal-insulator transition in 
the half-filled Hubbard model is one of the fundamen- 
tal topics in strongly correlated electron systems X The 
Hamiltonian of the Hubbard model is given as 



H 



fct 



-H.c.)- 



(i) 



where t is the hopping integral between neighboring sites, 
U is the on-site Coulomb interaction strength, and [i is 
the chemical potential maintaining the average particle 
density at half filling. This is achieved by setting /i = U/2 
due to the presence of the particle-hole symmetry. cJ CT 
(ci a ) denotes the creation (annihilation) operator of an 
electron with spin direction a (=t, i) on the i-th site and 
n-icr = c lcr c i<j is the number operator. Hereafter we use 
t = 1 as the unit of energy. 

On the square lattice, assuming the paramagnetic 
state, the critical Coulomb interaction for the metal- 
insulator transition (Z7 C ) is roughly equivalent to the band 
width, as has been estimated by various numerical meth- 
ods such as the cluster-extended dynamical mean-field 
theory (CDMFT)^ variational cluster approximation 
(VCA)^ and variational Monte Carlo method^ If the 
antiferromagnetism (AFM) is allowed, however, the sys- 
tem immediately goes into the antiferromagnetic insula- 
tor state by switching on U, i.e., Uaf — 0, due to the 
perfect nesting of the Fermi surface with the nesting vec- 
tor Q = (tt, 7r). 

On the honeycomb lattice, situation is different from 
that of the square lattice. The tight-binding band dis- 
persion exhibits a semimetallic (SM) behavior with the 
Dirac cones at the corners of the hexagonal Brillouin zone 
(K and K' points) and vanishing density of states at the 
Fermi level. The perfect nesting of the Fermi surface is 
absent due to the point-like Fermi surface and thus the 
critical Coulomb interactions for the semimetal-insulator 
transition (SMIT) and AFM transition can differ from 
those of the square lattice. 



So far, the quantum phase transitions in the honey- 
comb lattice have been investigated by various theoretical 
and numerical studies. A renormalization group anal- 
ysis in the large- N limit,— where N is the number of 
fermion flavor, predicted that, if the analysis holds in 
the physical case N = 2, the SM-AFM transition occurs 
at a finite U/t value. The Mott transition in a param- 
agnetic state has been studied by DMFT calculations^ 
where the SMIT is shown to occur at relatively large 
interaction strength (U c > 10). Since the honeycomb 
lattice is the two-dimensional lattice with the smallest 
coordination number (three neighbors), spatial correla- 
tions neglected in DMFT should be taken into account 
and thus CDMFT studies have also been performed. A 
finite-temperature CDMFT with continuous-time quan- 
tum Monte Carlo impurity solver anticipated U c = 3.3 
at zero temperature limits Another finite-temperature 
CDMFT with exact-diagonalization solver by LiebschiS 
showed that the critical interaction is about U c = 3 — 4, 
but the explicit value of U c was not found due to the diffi- 
culty from limitation of the energy resolution. Recently, 
He and LuAA showed that the SMIT occurs at U c = 
and AMF transition is of the first order with the lower 
and upper critical Coulomb interactions Uafi — 4.6 and 
Uaf2 = 4.85. Quantum Monte Carlo (QMC) simula- 
tions have also been performed by various authorsJ^r— 
In particular, the discovery of a spin-liquid phase be- 
tween U c — 3.5 and Uaf = 4.3 by Meng et aZ.jii has 
attracted much attention and stimulated researches on 
this issue. Recently, however, a large scale QMC up to 
2592 sites by Sorella et al. ruled out the existence of the 
spin liquid state and found a direct SM-AFM transition 
at Uaf = 3.76 ± 0.04 J£ All in all, critical behaviors near 
the SMIT and AFM of the Hubbard model on the hon- 
eycomb lattice are rather controversial than those on the 
square lattice. 

In this paper, motivated by such developments in the 
field, we study the AFM and SMIT in the Hubbard 
model on the honeycomb lattice by VCA with exact- 
diagonalization solverJ^r— We first calculate the AFM 



order parameter as a function of U with various clusters 
such as L c — 6, 8, 10 and (L c , L b ) = (6, 6) clusters, where 
L c and L\, are the number of correlated sites and bath 
sites, respectively (see FigfT]). Then, to study the SMIT, 
the single-particle excitation spectra and single-particle 
gap are calculated as a function of U. We thus show that 
U c is in fact vanishingly small and therefore the para- 
magnetic insulating phase appears in a wide interaction 
strength. We also calculate the spin correlation functions 
and find the development of localized moments, whereby 
the validity of effective spin models is discussed. We thus 
show that VCA is a suitable method to study the above 
issues not only because it can take into account the spa- 
tial correlations important in low-dimensional systems, 
but also because it can be performed at zero tempera- 
ture where the thermal fluctuations are absent and thus 
it can achieve a high energy-resolution essential for de- 
termining the SMIT. 
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FIG. 1: Clusters used in the present work, (a) L c — 8, 
(b) {L c ,L b ) = (6,6), (c) L c = 6, and (d) L c = 10. Filled 
(empty) circles represents the sites on the A (B)-sublattice of 
the honeycomb lattice and squares in (b) represents the bath 
sites. 



II. METHOD OF CALCULATION 

A. Variational cluster approximation 

We employ the VCAi£-i2 to study the SMIT and AFM 
of the Hubbard model on the honeycomb lattice at zero 
temperature. In order to study AFM, we introduce the 
staggered magnetic field 

Hy = ft/y^sign(z)o-rw, (2) 

where sign(i) = +1(— 1) for i e A(B)-sublattice on the 
honeycomb lattice. Moreover, for the (L c , Lb) = (6, 6) 
cluster, we introduce the bath hybridization 

Hv = V $>t ff a i(7 + H.c.) + e' a £ 4,0*, , (3) 

i.a i,a 

where aj CT (ai a ) denotes the creation (annihilation) op- 
erator of an electron with spin direction a on the i-th 
bath site and V is the hybridization between correlated 
sites and bath sites. The bath level is fixed at e' a = U/2 
so as to keep the average particle density of the corre- 
lated sites at half filling. Thus the fictitious magnetic 
field h! and hybridization strength V' are the variational 
parameters. The eigenvalue problems for the Hamilto- 
nian H 1 = H + + Hy* defined on the small clusters 
are solved and the single-particle Green's functions are 
calculated by the Lanczos exact-diagonalization method. 
Then the grand-potential functional Q,(h\ V) is calcu- 
lated from the Green's functions and the stationarity 
point (dfl/dh',dn/dV) = (0,0) is searched to find the 
physical solution. 

Note that there are other possible variational param- 
eters; e.g., enhancement of the hopping integrals be- 
tween the correlated sites (At) and the hopping inte- 
grals between the bath sites (tbath)- We have checked 



these variational parameters for some U values on the 
(L c , Lb) — (6, 6) cluster and found that the optimal val- 
ues are |Ai| < 0.15 and |ibath| < 0.15 and the optimal V 
values are affected only by ~ 0.01 \V'\ in the presence of 
these At or ibath , and moreover the change in the ground- 
state energy is negligible. Thus we omit the variations 
with respect to At and tbath- The similar discussion on 
the irrelevance of these variational parameters in the vari- 
ational processes has been made by Balzer et alA for the 
half-filled Hubbard model on the square lattice. 

B. Cluster perturbation theory 

We use the cluster perturbation theory (CPT)2ii~— to 
obtain the lattice Green's function 

Gg tt (k,u>) = 1- jr Gg PT (k,u;)e- k -( r *- r A (4) 

Here, a (= A, B) denotes the sublattice index, is the 
position of the i-th site in the cluster, and 

G CPT (k,Lj)=G'(uj)[I-V(k)G'(Lj)}- 1 (5) 

is the CPT Green's function^ - — where G'iw) is the 
Green's function matrix of the cluster and V(k) repre- 
sents the hopping matrix between the clusters. 

III. RESULTS OF CALCULATION 

A. Staggered magnetization 

Figure[2]shows the U dependence of the staggered mag- 
netization defined as 

m=~^T siga(i)a(n ia ), (6) 
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where (■ • •) denotes the ground-state expectation value 
calculated by use of the CPT Green's function^ A mean- 
field-theory (MFT) result reproduced from Ref. is also 
shown. We fine that Uaf strongly depends on the shape 
of the clusters, i.e., Uaf = 3.8, 1.5, 2.7, and 3.7 for L c = 
6, 8, 10, and (L c ,Lb) = (6,6), respectively. The order 
parameter increases continuously with increasing U from 
Uaf, and thus the transition is of the second-order. The 
characteristic linear increase of the order parameter m cx 
U - Uaf near U « U A f derived by the MFT-H is also 
observed in the VCA results. Focusing on the results for 
the L c = 6 and (L c , Lb) = (6, 6) clusters, we find that the 
bath degrees of freedom do not much affect the values of 
both Uaf and m. Our result Uaf = 3.7 obtained for the 
hexagonal cluster shows a remarkable agreement with the 
large-scale QMC simulation {U A f = 3.76 ± 0.04)±£ and 
reasonable agreement with the results of other previous 
studiesJaLr— On the other hand, Uaf calculated for the 
L c = 8 and 10 clusters are lower than that obtained by 
the MFT (Uaf = 3.2). Actually, being different from the 
L c = 8 and 10 clusters, the hexagonal clusters do not 
suffer from boundary effects and spatial homogeneity is 
kept as it should be in the thermodynamic limit. Thus 
we conclude that L c — 6 and (L c ,Lb) — (6,6) clusters 
are the relevant choice of the clusters and hereafter we 
show the results of this choice. 
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FIG. 2: (Color online) U dependence of the staggered mag- 
netization m calculated for various clusters (see also Fig. [TJ . 
The mean-field result (solid line) is reproduced from Ref. l2a . 



r 

M 
K 

r 
r 

M 
K 

r 




u=o 





U=2 




U=3 



t/=4 




U=5 



FIG. 3: Single-particle excitation spectra A(k, u) calculated 
for a paramagnetic state. Inset: definition of the F, M, and 
K points in the first Brillouin zone. 



paramagnetic state with the (L cl Lb) = (6,6) cluster are 
shown in Fig. [3J At U = 4, the spectra on the T-K and 
T-M lines become less dispersive, indicating the localized 
nature of electrons. On the contrary, on the K-M line 
(at the edge of the hexagonal Brillouin zone), the spec- 
tra remain to have quasiparticle-like sharp peaks even for 
U > 4. At U = 5, the upper and lower Hubbard bands 
appear clearly with a large intensity at uj ~ ±U. Al- 
though the single-particle gap exists clearly at U = 5, it 
is hard to judge when the gap opens by increasing U due 
to the artificial Lorentzian broadening 77. Thus we esti- 
mate the single-particle gap in a different way discussed 
in the next subsection. 



B. Single-particle excitation 

Next we calculate the single-particle excitation spectra 
defined as 

A(k,«) = --9f ]T G£ tt (k I w + tT7), (7) 

a=A,B 

where a small imaginary part r] of the complex frequency 
represents the Lorentzian broadening for the spectra. We 
use ry = 0.1. The calculated results from U = to 5 for a 



C. Single-particle gap 

The single-particle gap is defined as 

A=m + -/* _ , (8) 

where is the upper (lower) bound of the chemical 

potential. Since /i +< --^ can be evaluated by calculating 
the particle density by integrating the CPT Green's func- 
tion along the imaginary frequency axis^ A does not 
suffer from the artificial Lorentzian broadening. 

Another way to estimate the single-particle gap is to 
evaluate the intensity of the self-energy at the Fermi level. 
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FIG. 4: (Color online) (a) U dependence of the single-particle 
gap A calculated from Eq. ((HJ (triangles) and Eq. (|10p (in- 
verted triangles). (L c ,Lt) = (6,6) cluster is used for the 
calculation, (b) U dependence of the single-particle gap from 
Eq. (|10|l calculated for L c = 6 and (L c , Lb) = (6, 6) clusters. 




FIG. 5: (Color online) U dependence of the spin correla- 
tion function of the (L c , Lt) = (6, 6) cluster. Arrows indicate 
fiii, 5i3, Si4, and S12 (from top to bottom) for a 6-site Heisen- 
berg model with the nearest-neighbor exchange interaction. 
Inset: definition of the site index i. 



at infinitesimally small values of U, just like in the one- 
dimensional Hubbard model at half filling, i.e., U c = 0. 
To see the effect of bath sites, we show the results for the 
L c = 6 and (L c ,L b ) = (6,6) clusters in Fig. |U(b). For 
U > Uaf, results for both the AFM and paramagnetic 
states are shown. We find that the introduction of the 
bath sites significantly reduces the magnitude of the gap 
but it cannot close the gap, as we have seen in Fig. |U(a). 
For large U regime, the gap increases linearly with U, as 
it should be in the Mott insulator. 



The lattice self-energy £i att (k, lu) can be calculated from 
G^J t (k, ui) through the Dyson equation. The intensity of 
the self-energy at the K point on the Fermi level 

0"K = -w3£latt(K,Kj)| w ^ (9) 

is related to the single-particle gap in the following way. 
It is known that the self-energy has only the simple poles 
on the real frequency axis. Thus, if the single-particle 
gap exists, the self-energy near the Fermi-level has the 
form Ei a tt(K,cj) s» 2& Therefore the lattice Green's 
function at the K point near the Fermi level is approxi- 
mately given as Gi att (K, w) « \ {jj^— + s+^sj?) . Le -> 
the self-energy splits the spectrum above and below the 
Fermi level with the single- particle gap 

A ~ 2y^. (10) 

Figure 2] (a) shows the [/-dependence of the single- 
particle gap estimated from the chemical potential differ- 
ence Eq. © and intensity of the self-energy Eq. (fl"U| . Al- 
though the gap is extremely small (< O.OOlt for U < 0.2), 
the two results nicely coincide. We thus find that, even at 
U = 0.1, a small but finite single-particle gap exists. We 
may therefore conclude that the single-particle gap opens 



D. Spin correlation function 

Figure [5] shows the U dependence of the spin correla- 
tion function in the cluster 

S u - (*„|£fSf|*o>, (U) 

where |\&o) is the ground state of the (L c ,Lb) = (6,6) 
cluster with the optimal hybridization V' and Sf = 
y^ j(r <jni a /2 is the spin operator for the correlated site 
i. Positions of i-th sites are defined in the inset of Fig. 
[5] Calculations are done in the paramagnetic state. For 
comparison, the same quantities for a 6-site Heisenberg 
model with the nearest-neighbor exchange interaction are 
shown by arrows. 

The on-site correlation function Su represents the de- 
velopment of the local moment. It increases monotoni- 
cally and almost linearly with increasing U from Su = 
0.125 in the non-interacting limit toward Su — 0.25 in 
the localized spin limit. The second- and third- neighbor 
correlation functions S13 and 15141 take small values for 
U < 4. This indicates that, although the single-particle 
gap is finite, electrons still keep their itinerancy and spin 
correlations beyond the neighboring sites are not devel- 
oped, unlike in the Heisenberg model. Thus, to investi- 
gate the spin-liquid nature from an effective spin model, 
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the Heisenberg model with high t/U-order exchange in- 
teractions should be necessary. Antiferromagnetic cor- 
relations between the neighboring sites IS12I, as well as 
and I S'14 1 , start increasing from U ~ 4, where the 
electrons begin to localize and local moments are devel- 
oped. This is consistent with the emergence of the AFM 
at Uaf — 3.7 and the spectral signature of electron lo- 
calization (see Fig. 13]). In this interaction regime, the 
Heisenberg model with the nearest-neighbor exchange in- 
teraction may be appropriate for an effective low-energy 
model of the Hubbard model on the honeycomb lattice. 

IV. SUMMARY 

We have studied the semimetal-insulator transition 
(SMIT) and antiferromagnetism (AFM) in the half-filled 
Hubbard model on the honeycomb lattice within the vari- 
ational cluster approximation (VCA). The AFM transi- 
tion was found to be of the second order and the crit- 
ical Coulomb interaction for the AFM (Uaf) obtained 
in the hexagonal clusters showed a remarkable agree- 
ment with the recent large-scale quantum Monte Carlo 
simulation^ The single-particle gap has been calculated 



down to U — 0.1 and found to be finite. Thus we con- 
cluded that the critical Coulomb interaction for the SMIT 
is U c = 0. 

The extremely small gap obtained by our zero- 
temperature calculations in the small U regime suggests 
that the careful evaluation of the single-particle gap is 
required for the SMIT on the honeycomb lattice. In 
the paramagnetic insulating state, we have calculated 
the spin-correlation functions and found that the local 
moments are not developed. Thus we concluded that, 
the high t/ [/-order exchange interactions should be nec- 
essary to investigate the spin-liquid state from effective 
spin models. 
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